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ABSTRACT 

We develop a formalism for calculating probabilities for the outcomes of stellar dy¬ 
namical interactions, based on results from IV-body scattering experiments. We focus 
here on encounters involving up to six particles and calculate probabilities for direct 
stellar collisions; however our method is in principle valid for larger particle numbers. 
Our method relies on the binomial theorem, and is applicable to encounters involving 
any combination of particle radii. We further demonstrate that our base model is valid 
to within a few percent for any combination of particle masses, provided the minimum 
mass ratio is within a factor of a few from unity. This method is particularly suitable 
for models of collisional systems involving large numbers of stars, such as globular 
clusters, old open clusters and galactic nuclei, where small subsets of stars may reg¬ 
ularly have very close encounters, and the direct integration of all such encounters is 
computationally expensive. Variations of our method may also be used to treat other 
encounter outcomes, such as ejections and exchanges. 

Key words: gravitation - binaries (including multiple): close - globular clusters: 
general - stars: kinematics and dynamics - scattering - methods: analytical. 


1 INTRODUCTION 

The three-b ody pro b lem h as a long history extending all the 
way back to lNewtonl (1 16861 ). yet it has never been fully solved 
analytically. Indeed, for the majority of the relevant param¬ 
eter space and over suffi ciently long timescales, the evolu¬ 
tion is chaotic (IPoincarei 118921 ). More recently, a number of 
different appro aches have b een invoked to st udy the prob¬ 
lem (e.g. lHenonlll969l : IValtonen fc Karttunenll2006l ). and the 
most common approach today is to use computers to inte¬ 
grate the equations of motion direc tly (in some instances, us - 
ing the secular approximation, e.g. iNaoz fc Fabrvckvll2014l) . 
And yet, in modern IV-body simulations of star cluster evo¬ 
lution, it is the presence of binaries and especially higher- 
order multiples that poses the biggest computational chal¬ 
lenge. The time evolution of the orbits within multiple 
star systems, as well as the outcomes of chaotic close 
gravitational encounters involving multiples, requires very 
small time st eps, relative t o the crossing time of the entire 
cluster (e.g. iHurlev et al.l 120051 : I Geller, Hurley fe Mathieul 
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120131) . This is particularly problematic, since simulations 
have shown that encounters involving binaries and triples 
can be crucial to not only the overall cluster e volu- 
tion (e.g. iHutl 1 1983al : iHut. McMillan fc Romanil [i~992l) . but 


stragglers (Perets & Fabrvckv 

2009j: Leigh & Sills 

2011 

Geller. Hurley & Mathieu 201f 

1: INaoz & Fabryckv 

2014) 


more, observations have now revealed that higher-order mul¬ 
tiple star systems are p resent in young st ar clusters in non- 
negligible numbers (e.g. iLeigh et al.ll2013l ). 


One avenue into alleviating the computational chal¬ 
lenges of directly integrating these relatively low-A systems 
(where N is the number of particles involved in the interac¬ 
tion) is to assume that the final states specifying the out¬ 
comes of chaotic gravitational interactions are probabilistic, 
as opposed to deterministic, in so far as how they relate to 
the initial encounter conditions. For example, a strongly in¬ 
teracting three-body system ty pically b reaks up to produce 
a binary and an escaping star. IMonaghanl (I1976al ) showed 
that, by making the assumption that the probability of a 
configuration is proportional to the associated volume in 
phase-space, the statistical properties of the remaining bi- 
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nary and the escaper can be predicted (lMonaghan[ 197 6bli. 
This result rests heavily on Liouville’s Theorem (Liouvillc 
Il838h . which states that, in a gravitionally-interacting sys¬ 
tem of particles, the total volume in phase-space is conserved 
in time. _ 

In this, and our previous paper (Paper h lLeigh fe Geller] 
I2012IT we develop a probabilistic interpretation of the out¬ 
comes of stellar encounters involving modest numbers of par¬ 
ticles, with an over-arching goal of connecting the small- and 
large -N limits in bound configurations of finite-sized parti¬ 
cles evolving under the influence of Newtonian gravity. In 
Paper I, we performed numerical scattering experiments in¬ 
volving jsingle, binar y an d triple stars, using the FEWBODY 
code (IFregeau et alj|2004h . and showed that the probability 
of a collision occuring during a chaotic gravitational inter¬ 
action involving N identical particles is proportional to N 2 . 
Interestingly, this same IV-dependence is what is predicted 
by the mean-free path approximation in the limit of very 
large N. 

I 11 this paper, we delve further into the physical origin 
of this IV-dependence, which relates to the binomial theo¬ 
rem, and explain how to exploit it to predict collision prob¬ 
abilities over a large range in the relevant parameter space 
of initial conditions with minimal computational expense. 
Specifically, here we simulate encounters involving different 
types of particles (whereas in Paper I we limited our study 
to identical particles). We focus first on encounters involving 
particles with identical masses but different physical radii, 
and outcomes that result in a single direct collision between 
any two particles. We then go on to consider encounters 
involving particles with different masses. In Section [2] we 
present our method for calculating probabilities for the en¬ 
counter outcomes (focusing on collisions, but with a discus¬ 
sion of how to expand these results to different outcomes). 
Our results are presented in Section [3] which are then sum¬ 
marized and dicussed in Section [J 


2 METHOD 


path approximation. Specifically, in the limit of very large 
N, the mean free path approximation predicts that the rate 
of collisions between any two particles also scales as N 2 . The 
physical origin of this proportionality can be traced back to 
the binomial theoremlj 

Briefly, the binomial theorem states that the number 
of subsets composed of k particles that can be constructed 
from a larger sample composed of N particles is equal to: 


N\ _ N\ 
k ) ~ (N — k)\k\ 


( 1 ) 


Thus, in our specific example of collisions, we want the prob¬ 
ability of selecting any two particles at random from our 
larger sample of N identical particles. From the binomial 
theorem, the number of ways to select any two particles is: 

(N\_ N\ _ N{N - 1) 

(TV — 2)!2! 2 1 j 

This factor of N(N — 1) is proportional to the trend in colli¬ 
sion probability as a function of N particles that we found in 
Paper I (as stated above). Furthermore, N 2 ~ N(N — 1) for 
N 1, which is the scaling predicted by the mean free path 
approximation for collisions between any two particles. That 
is, the probability of any two identical particles colliding is 
proportional to Equation [2] In this way, the binomial theo¬ 
rem will form the backbone of any theory aiming to predict 
the relative probabilities of colliding two specific particles 
during a chaotic gravitational interaction. 

What if the particles are not identical? Specifically, let 
us consider the case of two different particles, types 1 and 2. 
Three different collision scenarios are now possible; a colli¬ 
sion between two particles of type 1, a collision between two 
particles of type 2 and a collision between a particle of type 
1 and a particle of type 2. Thus, the collision probability 
can be broken up into the sum of three different terms, each 
corresponding to the probability of a specific collision event 
occurring: 


Peon = P 11 + P 12 + P 22 , (3) 


In this section, we present our method for calculating col¬ 
lision probabilities for encounters involving different types 
of particles, as well as the numerical scattering experiments 
used to test our extrapolation technique. We will focus first 
on encounters where the particle masses are all identical, 
but two different particle radii are possible. However, we 
keep our approach as general as possible, since in principle 
it can be extended to any combination of particle numbers 
and types. Later we relax the constraint of equal particle 
masses and begin a comparison to encounters involving par¬ 
ticles of different masses. 


2.1 Model 

The key result from Paper I can be explained as follows. 
First, we showed that the probability of a collision occur¬ 
ing during a chaotic gravitational interaction involving N 
identical particles is proportional to N 2 . (Note that a pro¬ 
portionality of N(N—1) also fits the data; however we could 
not differentiate between these IV-dependences at a statisti¬ 
cally significant level.) We briefly touched upon a connection 
between this result and the large-IV limit via the mean free 


where P\ 1 and P 22 are the probabilities for a collision be¬ 
tween two particles of type 1 and two particles of type 2, re¬ 
spectively, and P \2 is the probability for a collision between 
a particle of type 1 and a particle of type 2. Now, the de¬ 
pendence of each of these probability terms depends on the 
number of particles of each type Ni, and this W-dependence 
comes from the binomial theorem. 

Let Ni and N 2 be the numbers of particles of type 1 
and type 2, respectively, that are involved in the dynamical 
encounter. Then, Equation [3] can be re-written as: 


Pcoll = Oil 





(4) 


where an, ai 2 and a 22 are constants that depend on the 


1 More technically, the binomial theorem gives the different pos¬ 
sible pairs of objects that could collide directly, and provides a 
good intuitive basis for building our model for the collision prob¬ 
ability. Hence, our model implicitly assumes that other effects 
relevant to deciding the outcome of a chaotic gravitational en¬ 
counter, and hence the collision probability, do not also depend 
on the number of interacting stars N. 
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total encounter energy E , the total angular momentum L 
and the properties of the particles (mass and radius). 

We further assume that the particles have identical 
masses but two different radii are possible, namely Ri and 
R 2 for particle types 1 and 2, respectively. We assume that 
the constant ai 2 can be found directly from the constants 
an and 022 : 

ail<Tll + 022(722 ,_ x 

012 = - ~ XT -. ( 5 ) 

where an and 022 are the collisional cross-sections for direct 
collisions between two particles of type 1 and two particles of 
type 2, respectively. The particles have identical masses, so 
the collisional cross-sections for these two collision scenarios 
are independent of mass and are set equal to the sum of their 
geometrical cross-sections, or an = 2-kR\ (for example). 
Hence, from Equation [5] the constant 012 is calculated as 
the mean of the constants an and 022 , weighted by the 
corresponding geometrical cross-sections. 

Equations [4] and [5] can be understood as follows. Dur¬ 
ing the evolution of a chaotic gravitational encounter, stars 
are continually undergoing close approaches with each other. 
The probability of a close approach corresponding to a sep¬ 
aration smaller than the sum of the particle radii occurring 
over the course of a crossing time is proportional to the num¬ 
ber of interacting stars. That is, over a given period of time, 
more clo se approaches occur if more particles are involved 
(IValtonen fe Karttunenll2006l j. However, if all particles are 
identical, then the probability of an event occurring (e.g. 
collision, ejection, etc.) in a crossing time is the same for 
all particles. Hence, identical particles have the same prob¬ 
ability of experiencing a close approach corresponding to a 
separation smaller than the sum of the particle radii (which 
we call a collision). If the particles are not identical, then 
the volume of phase-space corresponding to a collision de¬ 
pends on the types of the colliding particles. The constants 
an, a 12 and 0122 in Equation [S] correct for these different 
phase-space volumes for different collision scenarios, which 
we assume are (non-linearly) proportional to their respective 
collisional cross-sections. 

The key point is that, knowing a priori the constants 
an and 022 (he. from numerical scattering experiments), 
the collision probabilities for the remaining IV — 1 different 
interaction scenario^] can be calculated directly without the 
need for any additional numerical scattering experiments. 
As we will show in Section [3] for a given total encounter 
energy and angular momentum, a relatively small set of nu¬ 
merical scattering experiments can be extrapolated from to 
directly obtain the coefficients an and 022 for any particle 
radii. Additionally, numerical scattering experiments need 
only be run for the 2+2 or 1+3 case^f] since, using the bi¬ 
nomial theorem, the corresponding coefficients for higher-IV 

2 For a given N , there are N + 1 total encounter scenarios, as is 
evident upon looking at Figure [T] For the 2+3 case, this figure 
illustrates that there are six interaction scenarios in total (since 
five particles are involved in the interaction). The coefficients for 
the two cases involving all identical particles are obtained directly 
from numerical scattering experiments. Collision probabilities are 
then calculated for the four remaining interaction scenarios, with 
different numbers of stars of types 1 and 2, using these two coef¬ 
ficients. 

3 In principle, the 1+2 case can even be used, but here the nor- 


interactions (2+3, 3+3, etc.) can be calculated directly. This 
offers a strategy for using a minimal number of numerical 
scattering experiments to calculate collision probabilities for 
encounters involving any number of particles with any com¬ 
bination of particle radii. 

2.2 Numerical scattering experiments 

In this section, we present the numerical scattering experi¬ 
ments used to test our strategy for calculating collision prob¬ 
abilities presented in the previous section. 

We calculate the outcomes of a series of binary-binary, 
single-triple, binary-triple, and triple-triple encounters using 
the FEWBODY numerical scattering codcQ. The code integrates 
the usual N-body equations in configuration- (i.e., position- 
) space in order to advance the system forward in time. 
For details concerning the adaptive integration, classifica- 
tio n tec hni ques, etc, used by FEWBODY, we refer the reader to 
iFregeau et al.l (|2004| ~). 

In Paper I, we adapted the FEWBODY code to handle en¬ 
counters involving not only single a nd b inar y stars , b ut a lso 
triples[f] We use the same criteria as lFregeau et al.l (120041 ') to 
decide when a given encounter is complete, defined as the 
point at which the separately bound hierarchies that make 
up the system are no longer interacting with each other or 
evolving internally. 

To perform physical collisions between stars, FEWBODY 
uses the “sticky-star” approximation. This treats stars as 
rigid spheres with radii equal to their stellar radii. A physi¬ 
cal collision is assumed to occur when the radii of the stars 
overlap. When this happens, the stars are merged assuming 
conservation of linear momentum and no mass loss. This 
does not account for tidal effects, whi ch could signific antly 
increase the collision probability fe.g. lMcMillanlfl986lj . but 
are beyond the scope of this work. For this reason, we con¬ 
sider the collision probabilities presented in this paper to be 
lower limits for the true values. 

Previous scattering experiments have shown that the 
total energy and angular momentum are critical param- 
eters in deci ding the outco mes of 1+2 interactions (e.g. 
IValtonen fe Karttunenll2006r ). We confirmed this for higher- 
N interactions (up to N = 6) in Paper I. Therefore, we 
will fix the total energy and angular momentum when com¬ 
paring between encounters involving different numbers of 
objects. By fixing these quantities, we aim to remove the 
dependences of the encounter outcomes on energy and an¬ 
gular momentum (when comparing between encounters in¬ 
volving different numbers of particles, or different combina¬ 
tions of particle types/radii), thereby normalizing the com¬ 
parisons to reveal the dependences of the collision probabil¬ 
ity on other parameters, in particular the particle radii. We 

malization in tota l energy and angular m omentum is slightly more 
complicated (see iLeigh &; Gelled 120121 for more details), since 
there is only one binary orbit going into the interaction. 

4 for the source code, see http://fewbody.sourceforge.net 

5 Specifically we created additional subroutines to simulate 1+3 
and 3+3 encounters; codes to simulate encounters between bina¬ 
ries and singles only, as well as a 2+3 encounter code, were pre¬ 
viously available in the FEWBODY package. The authors are happy 
to provide these additional subroutines to users of FEWBODY upon 
request. 
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consider one specific combination of the total energy and 
angular momentum (see below), for which both the total 
energy (range) and angular momentum are always chosen 
to be the same to within a factor of ~ 2 for every type of 
encounter (see Paper I for more details). 

To this end, all orbits are circular, and have semi-major 
axes of either ao = 0.1 AU or ai = 1 AU. Hence, for a 2+3 
encounter for example, we adopt ao = 0.1 AU for the or¬ 
bital separation of the inner orbit of the triple, and ai = 1 
AU for both the outer orbit of the triple and the interloping 
binary. For a 1+3 encounter, we similarly adopt ao = 0.1 
AU and ai = 1 AU for the inner and outer orbits, respec¬ 
tively, of the triple, while for a 2+2 encounter we adopt ao = 
0.1 AU and ai = 1 AU for the two binaries. These initial 
semi-major axes are summarized in Table [T] A semi-colon is 
used to separate different objects, whereas a comma is used 
to separate the orbits within triples. Parantheses are used 
to enclose the semi-major axes of triples, with the smaller 
of the two semi-major axes always corresponding to the in¬ 
ner binary. Note that we have checked that our results are 
approximately insensitive to our choice of semi-major axes, 
for orbital separation ratios ao/ai < 100. For a given rel¬ 
ative velocity at infinity and impact parameter, the choice 
of semi-major axes sets the encounter energy and angular 
momentum. Thus, while the exact collision probability does 
depend on the choice of semi-major axes, our method does 
not. This choice of semi-major axes, along with our cho¬ 
sen impact parameter (6 = 0; see below), serves to minimize 
the computational expense of our scattering experiments, by 
minimizing the total angular momentum. Regardless, these 
semi-major axes are relatively characteristic of the multiple 
star systems we expect to be present in globular clusters. At 
least in the cluster core, only relatively compact multiples 
should survive for several Gyr in such dense stellar environ¬ 
ments. This is especially true for triples (and higher-order 
multiples), given that they are only stable against disso¬ 
ciation if the ratio between their i nner and outer orbital 
separations is sufficiently small fe.g. iMardlindEoQll 'l. 

For each Run, we populate a grid of scattering exper¬ 
iments varying only the relative velocity at infinity i>i n f. 
Specifically, we select values for the relative velocity at in¬ 
finity from 0 to 1.1 in equally-spaced intervals of 0.004, in 
units of the critical velocity u cr it (defined as the relative ve¬ 
locity that gives a total energy of zero for the encounter). 
Hence, our chosen range of relative velocities at infinity in¬ 
cludes the case of zero total energy. In the case Vi n f > v cr it, 
the encounter is prompt and at least one particle is almost 
immediately lost from the system. For Vi n f < v cr it, all par¬ 
ticles may be temporarily bound and participate in a pro¬ 
longed resonant encounter. We fix the impact parameter at 
6 = 0, independent of Vi n {. For our assumed circular orbits, 
this minimizes the total angular momentum (and hence re¬ 
duces the total computational cost of our simulations) while 
preserving our normalization that the total energy and an¬ 
gular momentum must be roughly the same independent of 
the number of interacting particles. Finally, at each point in 
this grid, we perform multiple scattering experiments ran¬ 
domizing all angles in the encounter, including the angles 
at impact between the orbital planes of any binaries and 
triples involved in the encounters. In total, we perform 24750 
simulations for every Run (with the exception of Figure [4] 
for which we are limited to using only 4500 simulations per 


data point). We chose this number of simulations to balance 
between statistical significance and computational expense 
(i.e. simulation run-time). Of course, running more simula¬ 
tions would reduce the uncertainties on the calculated colli¬ 
sion probabilities and yield more precise results, which may 
be desirable for certain applications. 

The collision probability is calculated from the out¬ 
put of our numerical scattering experiments for a given en¬ 
counter type and a given Run as: 


Pcoii — 


Well 

Wot ’ 


(6) 


where Woli is the total number of encounters that result in 
a direct physical collision, and Wot is the total number of 
encounters performed (throughout this paper, Wot = 24750, 
except for Figure [4] as stated above). The uncertainty for 
the collision probability is calculated using Poisson statistics 
according to: 


APcoll = ((^jf ) 2 + APoVint) 172 , (7) 

where APcoii^nt is an additional term to account for any 
intrinsic dispersion in our calculated collision probabilities 
(see Section [3] for our assumed values of AP co n,mt). This 
additional source of dispersion comes from a combination 
of numerical inaccuracies, our normalization in energy and 
angular momentum, and subtle inadequacies of our model 
in capturing all of the relevant underlying physics (e.g. the 
issue of saturation which we discuss in more detail in Sec¬ 
tion O 0 Technically, this uncertainty should include scat¬ 
tering experiments that result in unresolved outcomes dHutl 
Il983bh . However, we find that the number of unresolved 
outcomes is sufficiently small that it does not significantly 
contribute to AP co n, and we do not include it in its calcu¬ 
lation. 


3 RESULTS 

In this section, we first present some example applications 
of our method. Using these, we go on to describe how to 
calculate collision probabilities for any number of particles 
and any combination of particle radii, as well as for a limited 
range of particle masses. 


3.1 Different particle radii 

First, consider the case where the particles all have the same 
mass, but two different particle radii are possible. We assume 
Pi = 0.1 Rq and P 2 = 5 Rq, and mi = m 2 = 1 Mg. Fig¬ 
ure □ shows the results of a series of numerical scattering 
simulations of binary-triple encounters (i.e. N = 5). Here, 
we vary the number of large particles from 0 to 5, as shown 
on the z-axis. The constants <rn and 022 in Equation [4] are 
calculated from our numerical scattering experiments using 
the W = 5 and W = 5 cases, where all particles are either 


6 Formally, the intrinsic dispersion term AP co n ; nt should be 
found by forcing the reduced chi-square of a fit to be unity. How¬ 
ever, in our case, we have only three data points (i.e. IV = 4, 5 and 
6) to work with in calculating the intrinsic dispersion term. Thus, 
necessarily, our estimates for AP co n ; nt are only approximate. 
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Table 1 . Initial semi-major axes of all binaries and triples 


Encounter Type 

Semi-major axes 
(in AU) 

2+2 

0.1; 1.0 

1+3 

(0.1, 1.0) 

2+3 

1.0; (0.1, 1.0) 

3+3 

(0.1, 1.0); (0.1, 1.0) 


of type 1 (an) or type 2 (022)- The constant ai2 is then 
calculated from Equation 0 The open points show the re¬ 
sults of these calculations, and are compared directly to our 
numerical scattering experiments via the filled points. Error 
bars are calculated using Equation [7] with AP co n.i n t = 0.05, 
which is chosen to ensure that the residuals calculated in 
Figure 0 are all within one standard deviation of zero. Our 
calculated predictions agree very well with the results of our 
numerical scattering experiments, and are within one stan¬ 
dard deviation of each other in all cases. The agreement is 
slightly worse for the N 2 = 1 and N2 = 2 cases, however, 
for which our calculated collision probabilities deviate from 
our numerical scattering experiments by about one standard 
deviation (see the bottom panel in Figure0. The reason for 
this poor agreement is not understood. This may be high¬ 
lighting that we are missing some higher-order term(s) in 
Equation 31 possibly related to our calculation of the coef¬ 
ficient (J12. Regardless, the calculated collision probabilities 
(via Equation 0 always agree with the results of our nu¬ 
merical scattering experiments to within a couple percent 
or better. 

In Figure [2] we show the results of a series of 2+2 (filled 
circles), 1+3 (open circles), 2+3 (filled triangles) and 3+3 
(open stars) encounters, this time involving all identical par¬ 
ticles. The encounter parameters are the same as in Figure [Tj 
and Table 0 with the sole exception that the particle radii 
are all identical. For each encounter type, we show the col¬ 
lision probability as a function of the ratio R/ao, where R 
is the particle radius and ao is the initial semi-major axis of 
the most compact orbit going into the encounter (i.e. a 0 = 
0.1 AU). Hence, R/ao = 1 corresponds to a contact config¬ 
uration for the most compact orbit initially. 

The key point is that an analogous version of Figure 0 
can be made for any combination of particle radii (i.e. any 
Ri and ifo). In principle, this can be done by extracting the 
coefficients <rn and 022 directly from Figure [2] and fitting a 
curve to these data to derive the collision probability at any 
value of R (and N). The temptation is, for every N, to fit 
the data with a function of the form: 

Fcoii = SR' 1 (8) 

However, as is clear from Figure 0 a plot of R versus P C oii 
is clearly sub-linear. Even in log-log space, significant curva¬ 
ture is apparent. This is not unexpected, since in the limits 
R —> 0 and R —> ao (where ao is the semi-major axis of the 
most compact orbit going into the encounter) we expect, 
respectively, P co ii —> 0 and P co ii —> 1, independent of the 
number of interacting particles N. Thus, the coefficients 5 
and 7 in Equation [8] are not constants, but are instead func¬ 
tions of the particle radius R. (In practice, one could fit a 


function similar to Equation [8] to a subset of data around 
the desired value of R/ao , in order to derive the collision 
probability at a point not plotted in Figure 0. 

Since here, by construction, the total encounter energy 
and angular momentum are roughly the same independent 
of N, then, via the binomial theorem, at a given value of 
the particle radius R, S and 7 should be the same for all 
N. That is, since the particles are all identical and the total 
encounter energy and angular momentum are fixed, the col¬ 
lision probability should vary only with the number of inter¬ 
acting particles N. To illustrate this point, we calculate the 
quantity -Pcoii/(^f) = an for encounters involving all identi¬ 
cal particles, for all values for the particles radius, R, shown 
in Figure 0 with the expectation that the ratio P c oii/(/f) 
will be the same for all N (to within the uncertainties). The 
results are shown in Table0and Figure0 Uncertainties are 
calculated using Equation0with APcdi^nt = 0.01, which en¬ 
sures that our simulated collision probabilities for the 2+2 
and 1+3 cases are always within one standard deviation of 
each other. 

The quantity P C oii/(/f) is approximately the same (e.g., 
within typically less than one standard deviation) for all 
IV at a given particle radius R, until reaching a “satura¬ 
tion” limit at log(P/ao) <) —1.5. This suggests that, at least 
for a given particle radius R, the collision probabilities for 
2+3 and 3+3 encounters can be calculated directly from 
the collision probability for the 2+2 (or 1+3) case, with¬ 
out having to run any 2+3 and 3+3 numerical scattering 
experiments. This offers the potential to calculate collision 
probabilities for high-A interactions with minimal compu¬ 
tational expense, which increases steeply with increasing N 
(as we also point out in Paper I). Equivalently, for fixed 
computational expense a larger number of simulations can 
be performed for lower N interactions, which subsequently 
increases the statistical significance and accuracy of the re¬ 
sults. We caution that the simulated data begin to deviate 
from the trend that Pcoii/^) is roughly independent of N 
near the saturation limit P co n —> 1. As shown in Figure 0 
this begins around log(P/ao) ~ —1.5, and should be treated 
with caution. More work is needed to verify that this exact 
saturation limit is ubiquitous for all total encounter ener¬ 
gies and angular momenta. In general, however, adopting 
larger initial semi-major axes a 0 pushes the saturation limit 
to larger values of the particle radius R. 


3.2 Different particle masses 

Next, we consider encounters involving particles having dif¬ 
ferent radii and different masses. The key point we wish 
to illustrate in this section is that, provided the range of 
energies and angular momenta are fixed and independent 
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Table 2. The quantity P co \\/ (^f) as a function of the particle radius R 


Particle Radius 
(in R 0 ) 

(2+2; N= 4) 

Pcon/O 

(1+3; iV=4) (2+3; N= 5) 

(3+3; N=6 ) 

0.005 

0.0042 

± 

0.0017 

0.0059 ± 0.0017 

0.0052 ± 0.0010 

0.0060 

± 0.0007 

0.01 

0.0069 

± 

0.0017 

0.0085 ± 0.0017 

0.0087 ± 0.0010 

0.0089 

± 0.0007 

0.05 

0.0203 

± 

0.0017 

0.0231 ± 0.0017 

0.0219 ± 0.0010 

0.0213 

± 0.0007 

0.1 

0.0304 

± 

0.0017 

0.0333 ± 0.0017 

0.0307 ± 0.0011 

0.0278 

± 0.0007 

0.5 

0.0618 

± 

0.0018 

0.0633 ± 0.0018 

0.0524 ± 0.0011 

0.0434 

± 0.0007 

1.0 

0.0773 

± 

0.0018 

0.0770 ± 0.0018 

0.0610 ± 0.0011 

0.0481 

± 0.0008 

5.0 

0.1255 

± 

0.0019 

0.1252 ± 0.0019 

0.0785 ± 0.0011 

0.0564 

± 0.0008 

10.0 

0.1581 

± 

0.0020 

0.1567 ± 0.0020 

0.0896 ± 0.0012 

0.0667 

± 0.0008 

20.0 

0.1667 

± 

0.0020 

0.1667 ± 0.0020 

0.1000 ± 0.0012 

0.0667 

± 0.0008 
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Figure 1. The probability of a collision occurring during a 
binary-triple encounter, shown in the top panel. Two particle 
radii are considered, namely R\ = 0.1 R© and R 2 = 5 R©. The 
number of large particles is shown on the z-axis, and the colli¬ 
sion probability is shown on the y-axis. All particle masses are 
identical, and the range of total encounter energies and angular 
momenta over which we integrate are identical for all cases (i.e. 
combinations of small and large particles). We adopt ao =0.1 
AU for the orbital separation of the inner orbit of the triple, and 
a± = 1 AU for both the outer orbit of the triple and the inter¬ 
loping binary. All orbits have eccentricities of zero, initially. The 
open triangles show our predictions for the collision probabilities, 
calculated from Equations [4] and [5] using the results of our nu¬ 
merical scattering experiments for the N± = 0 and N 2 = 5 cases 
(i.e. all small and all large particles). The filled triangles show 
the same collision probabilities, but calculated directly from the 
results of the corresponding numerical scattering experiments. Er¬ 
ror bars are calculated using Poission statistics (filled triangles) 
and simple error propagation (open triangles), beginning with 
Equation [7] and assuming AP co n j nt = 0.05. The bottom panel 
shows the residuals between the collision probabilities calculated 
directly from our numerical scattering experiments (filled trian¬ 
gles) and our model predictions (open triangles). The dashed line 
is shown for reference, and represents perfect agreement between 
our model and the data. 


Figure 2. The logarithm of the probability of a collision occur¬ 
ring during encounters involving all identical particles is shown 
as a function of the logarithm of the ratio R/ao , where R is the 
particle radius in AU and ao = 0.1 AU is the initial semi-major 
axis of the most compact orbit going into the encounter. Hence, 
we adopt ao = 0.1 AU for the orbital separations of the inner 
orbits of any triples, and ai = 1 AU for the outer orbits of any 
triples. For the 2+3 case, the initial semi-major axis of the bi¬ 
nary is ai = 1 AU, whereas for the 2+2 case the binaries have 
different initial semi-major axes of ao =0.1 AU and a± = 1 AU. 
Thus, ao is a constant for all encounters, independent of the total 
number of particles or the particle radii. The filled circles, open 
circles, filled triangles and open stars correspond to 2+2, 1+3, 
2+3 and 3+3 encounters, respectively. Error bars are calculated 
using Equation 0 with AP co n j nt = 0.01. 

of N , Equation [3] is approximately valid and describes the 
collision probabilities to within a few percent, independent 
of the distribution of particle masses, provided they are all 
within a factor < 2. Figure [4] shows two such examples. The 
first is shown by the open triangles, for which the particle 
masses are fixed and independent of our choices for the par¬ 
ticle radii. We set the particle masses to 0.7, 0.8, 0.9, 1.0 
and 1.0 Mq, with particle radii of either Ri = 0.1 R© or 
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Figure 3. The logarithm of the coefficient an given in Table [2] 
calculated as P C oll/(^)i a function of the logarithm of the 
ratio R/ao, where R is the particle radius in AU and ao = 0.1 
AU is the initial semi-major axis of the most compact orbit going 
into the encounter, for encounters involving all identical particles. 
We assume the same encounter parameters and symbols as in 
Figure [2] Error bars are calculated as in Figures [l] and [2] 

R 2 = 5 Rq assigned randomly (i.e. independent of particle 
mass). The initial configuration (i.e. which objects form the 
inner binary of the triple, etc.) is also chosen at random. 
The second example is shown by the open squares. Here, we 
vary the particle masses along with the particle radii, such 
that particles of type 1 have Ri =0.1 Rq and mi =0.5 Mg, 
while particles of type 2 have R2 = 5 Rg and m2 = 1 Mg. 
Again, the initial configuration is then chosen at random. 
Figure [4] shows that in both of these examples, Equation [3] 
accurately describes the collision probabilities to within a 
few percent. We have checked that this is also the case for 
mass ratios ~ 0.05, however a more thorough exploration of 
the relevant parameter space is required to accurately de¬ 
termine to what mass ratios this approximation holds. Our 
results support the conclusion that it is valid for mass ratios 
> 0.5, and possibly much lower. 

3.3 Putting it all together 

In this section, we combine the results of the previous two 
sections to construct a simple procedure for calculating col¬ 
lision probabilities in realistic star clusters, for any combina¬ 
tion of particle radii and masses. Thus far, we have focussed 
on encounters involving only two different types of particles. 
The method presented below, however, is more general. 

To summarize, our method is as follows: 

• Decide upon the encounter type and, in particular, the 
total number of particles involved in the interaction. As a 
guide, the observed fractions of binary and multiple star 
systems can be used to calculate analytic estimates for the 


Figure 4. The probability of a collision occurring during a 
binary-triple encounter. We assume the same encounter param¬ 
eters and particle radii as in Figure [T] This time, however, the 
particle masses can be different. The open and filled triangles 
show the same collision probabilities as in Figure [l] and assume 
all equal mass particles (i.e. 1 Mg). That is, the open triangles 
are calculated from Equations [4] and [5] using the results of our nu¬ 
merical scattering experiments for the A r i = 0 and IV 2 = 5 cases 
only, and the filled triangles show the same collision probabilities 
but calculated directly from the results of the corresponding nu¬ 
merical scattering experiments. The open circles show the same 
encounter parameters as in Figure |T] and are calculated directly 
from our numerical scattering experiments, but assuming parti¬ 
cle masses of 0.7, 0.8, 0.9, 1.0 and 1.0 Mg for every simulation. 
The open squares again show the same encounter parameters and 
are calculated directly from our numerical scattering experiments, 
but adopting radii and masses of Ri =0.1 Rq and m 1 = 0.5 Mg 
for particles of type 1, while particles of type 2 have radii R 2 = 
5 Rg and masses m 2 = 1 Mg. Error bars are calculated as in 
Figure [Til i.e. using Equation [7] with AP co u j nt = 0.05). 


rate s of the different encounter types (e.g. 2+2, 2+3, etc., 
see iLeonar dl ll989l : iLeieh fe Sillsll201 ll) . 

• Determine the initial total encounter energy and an¬ 
gular momentum from the initial orbital eccentricities and 
semi-major axes, along with the initial relative velocity at 
infinity and impact parameter. Note that the method is valid 
both for a particular combination of energy and angular 
momentum, as well as consistent ranges in these quanti¬ 
ties. Depending on how much information is available, one 
might estimate the relative velocities at infinity and the 
impact parameter in a star cluster from the obse rved sur¬ 
face brightness and velocity d ispersion profiles fe.g. iLeonardl 
Il989l : iBahramian et al.l l2013). Likewise, binary orbital pa¬ 
rameters and masses can be drawn from the appropriate 
observed (or assumed) distribution functions. 

• Set the particle types, in particular the particle radii 
Ri, along with the particle masses. Note that, as shown in 
Section El our method is valid to within a few percent for 
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any combination of particle masses (to within a factor of a 
few)Q 

• Given the particle radii, obtain the coefficients an, cry, 
etc. directly from numerical scattering experiments, using 
the parameters determined above. With these, calculate the 
coefficient (s) ay using the generalized equation: 


QiiOn + ay cry 
+ Ojj 


(9) 


where ay = 2nRf. As shown in Figure [5] nearly all possible 
values of the coefficients an, ay, etc. can be covered using 
a limited number of numerical scattering experiments by 
extrapolating between data points. Additionally, numerical 
scattering experiments need only be run for the 2+2 or 1+3 
cases since, using the binomial theorem, the corresponding 
coefficients for higher-TV interactions (2+3, 3+3, etc.) can be 
calculated directly (e.g., see end of Section T3. II) . This signifi¬ 
cantly reduces the total computational expense for higher -N 
interactions, since the integration times for the simulations 
increase dramatically with increasing particle number. 

• Calculate the desired collision probabilities using the 
generalized formula: 


Nj 

Pcoll — ^ ' O; 
i=l 


, Nj Nj 

+ ^E E 


i=i 


' Ni 

-m | 1 


( 10 ) 


As a proof of concept of the above generalized method- 
olgy, consider 2+3 encounters involving three different par¬ 
ticle types characterized by their different radii, or Ri = 
0.05 Rq, f?2 = 0.5 Rq and R 3 = 5.0 Rq. All particles have 
the same mass of 1 Mg, and the semi-major axes and ec¬ 
centricities are chosen to be the same as used in Figure [2] 
and Table [2] as is the range of relative velocities at infinity 
(including a small fraction of encounters with positive to¬ 
tal energies). As before, the impact parameter is set equal 
to zero for all simulations. We consider five different com¬ 
binations of particle types, and run 24750 simulations for 
each such combination. For comparison, we calculate the 
corresponding collision probabilities (along with their un¬ 
certainties erP co n) using Equations [177] and [7] as well as the 
coefficients an and ay obtained from Tablc[2] Based on Fig¬ 
ure 03 R3 = 5 R© is beyond the saturation limit for our 
model. This will introduce some small amount of additional 
scatter into our calculations for the collision probabilities. 
To account for this additional source of uncertainty, we cal¬ 
culated the mean value of an at R = 5 R© in Figure [3] for 
all values of the particle number N (i.e. N = 4, 5 and 6), 
and added the resulting value to its uncertainty (in quadra¬ 
ture) using Equation [7] with AP co u : i nt = 0.03. This serves to 
illustrate the importance of accurately estimating the coef¬ 
ficients an, ay, etc. when calculating collision probabilities. 

The results of this example are shown in Tabic [3] The 


7 In a real star cluster, the distribution of stellar masses and 
radii is continuous. Hence, some binning by stellar mass (which is 
tightly correlated with stellar radius on the main-sequence) is de¬ 
sirable, wherever this assumption is reasonable. For example, bin¬ 
ning is a reasonable assumption for interactions involving white 
dwarfs, bright red giants and/or main-sequence stars close to the 
turn-off. Here, the objects all have very similar radii and masses 
within their individual stellar evolution-based classifications, but 
very different radii between these classifications. 


agreement between our simulated and calculated collision 
probabilities is quite good in all cases, which generally agree 
to within a few percent. 


4 SUMMARY AND DISCUSSION 

In this paper, we present a formalism for calculating col¬ 
lision probabilities during chaotic gravitational encounters 
with finite-sized particles. The interactions involve any num¬ 
ber of particles and any combination of particle radii, using 
a minimal number of numerical scattering experiments. Our 
method relies on the binomial theorem, and can also be ex¬ 
tended to encounters involving particles of different masses 
(see below). 

Our model is applicable for a constant value or range in 
the total encounter energy and angular momentum. Hence, 
our formalism for the collision probability can be extended 
to include any combination of total encounter energy and 
angular momentum. To do this, plots analogous to Figure [2] 
are generated for different total encounter energies at fixed 
angular momentum (and vice versa). From here, the depen¬ 
dence of the fitting parameters (see Equation [8j) on the total 
encounter energy and angular momentum is derived. In gen¬ 
eral, the collision probability decreases with increasing total 
encounter energy (i.e. becomes less negative and thus closer 
to zero), and with increasing total angular momentum. De¬ 
pending on the difficulties associated with any curvature in 
log-log space (which prohibits fitting a constant power-law 
to the data), an analytic formula for the collision probability 
can be derived by running additional numerical scattering 
experiments to cover the full range in encounter energies 
and angular momenta. The resultant formula would thus 
cover the entire range of the relevant parameter space (to¬ 
tal encounter energy and angular momentum, distribution of 
particle masses and radii, particle number, etc., in, for exam¬ 
ple, open, globular or nuclear star cluster clusters). Alterna¬ 
tively, our method can be used to create a library of collision 
probabilities, again covering the entire range of the relevant 
parameter space while minimizing the total number of nu¬ 
merical scattering experiments that need to be performed 
(by extrapolating between data points, as in Figure [2]). This 
will be the focus of future work. 

Our formalism is valid as well for different particle 
masses, provided they are all within a factor of ~ a few. 
This is because changing the particle masses also changes 
the total encounter energy and angular momentum. Thus, 
by fixing the total energy and angular momentum for all 
encounters, the dependence of the collision probability is re¬ 
moved, at least to first-order. There are second-order effects 
related to the lifetime of the system, causing slightly higher 
collision probabilities for encounters with large particle mass 
ratios. However, the deviation is by at most a few percent 
from the model presented in this paper, which assumes all 
equal mass particles. 

Our method is ideal for application in old open and 
globular clusters, since it is applicable over a wide range 
of particle radii but a relatively narrow range of particle 
masses. Suitable interactions occur in GCs when single, bi¬ 
nary and triple stars undergo encounters, since the objects 
involved include main-sequence stars (close to the turn-off), 
red or asymptotic giant branch stars, white dwarfs and neu- 
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Table 3. The collision probability during 2+3 encounters for various combinations of particle types 1, 2 and 3 


Number of Each Type of Particle 

Collision Probability 

IVi 

n 2 

n 3 

Simulated 

Calculated 

Ri = 0.05 R 0 

R 2 = 0.5 Rq 

i?3 = 5.0 Rq 

Peon -k- ^coii 

-^coll i ^coll 

4 

1 

0 

0.3324 ± 0.0042 

0.3398 ± 0.0074 

4 

0 

1 

0.5264 ± 0.0057 

0.4454 ± 0.1397 

3 

2 

0 

0.4084 ± 0.0048 

0.4307 ± 0.0073 

3 

1 

1 

0.5670 ± 0.0060 

0.5357 ± 0.1103 

1 

1 

3 

0.7221 ± 0.0071 

0.7578 ± 0.1807 


tron stars. All of these objects have masses that are typically 
the same to within a factor < 2, but radii differing by several 
orders of magnitude. Thus, the geometric cross-section dom¬ 
inates the total collisional cross-section for these objects, as 
opposed to the gravitationally-focused cross-section. For ex¬ 
ample, consider an average encounter in a typical GC, in¬ 
volving mostly white dwarfs. If even a single main-sequence 
star is included in the interaction, then the collision prob¬ 
ability can increase by nearly an order of magnitude (this 
is the case for, e.g. the encounter parameters assumed in 
Figure [5] for which the collision probability begins to satu¬ 
rate at the particle radius corresponding to that of a typical 
main-sequence star in an old globular cluster) [f] Consider 
a 2+2 interaction with the same encounter parameters as 
adopted in Figure [2] and Table [2] If the interaction involves 
only white dwarfs, then P co ii ~ 0.04. If one of the white 
dwarfs is replaced by a single (near turn-off) main-sequence 
star, then P co ii ~ 0.25. Conversely, consider the same en¬ 
counter but involving only main-sequence stars with masses 
close to the turn-off mass. The collision probability for this 
encounter is P co ii ~ 0.46, roughly an order of magnitude 
larger than calculated for the same encounter but involv¬ 
ing only white dwarfs. Here, if a single main-sequence star 
is replaced by a white dwarf, then the collision probability 
changes to P co ii ~ 0.31. 
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